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It is generally assumed 1 ' 2 ' 3 that solid hydrogen will 
transform into a metallic alkali-like crystal at suf- 
ficiently high pressure. However, some theoretical 
models 4,5 have also suggested that compressed hy- 
drogen may form an unusual two-component (pro- 
tons and electrons) metallic fluid at low temperature, 
or possibly even a zero-temperature liquid ground 
state. The existence of these new states of mat- 
ter is conditional on the presence of a maximum in 
the melting temperature versus pressure curve (the 
'melt line'). Previous measurements 8 ' 7 ' 8 of the hy- 
drogen melt line up to pressures of 44 GPa have led 
to controversial conclusions regarding the existence of 
this maximum. Here we report ab initio calculations 
that establish the melt line up to 200 GPa. We pre- 
dict that subtle changes in the intermolecular interac- 
tions lead to a decline of the melt line above 90 GPa. 
The implication is that as solid molecular hydrogen 
is compressed, it transforms into a low-temperature 
quantum fluid before becoming a monatomic crys- 
tal. The emerging low-temperature phase diagram 
of hydrogen and its isotopes bears analogies with 
the familiar phases of 3 He and 4 He, the only known 
zero-temperature liquids, but the long-range Coulom- 
bic interactions and the large component mass ratio 
present in hydrogen would ensure dramatically dif- 
ferent properties—ii— . 

The possible existence of low-temperature liquid 
phases of compressed hydrogen has been rationalized 
with arguments based on the nature of effective pair in- 
teractions and of the quantum dynamics at high density, 
resulting in proton-proton correlations insufficient for the 
stabilization of a crystalline phased. But so far there 
has been no conclusive evidence establishing whether 
hydrogen metallizes at low temperature as a solid (the 
more widely accepted view to date) or as a liquid. Mea- 
surements and theoretical predictions of the near-ground 
state high-pressure phases^ of hydrogen have proven to 
be difficult because of the light atomic mass, signifi- 
cant quantum effects and strong electron-ion interac- 
tions. In this regard, the finite temperature liquid-solid 
phase boundary predicted here is especially valuable for 
understanding the manner in which hydrogen metallizes. 

The appearance of a maximum melting temperature in 
hydrogen is in itself a manifestation of an unusual phys- 
ical phenomenon. The few systems with a negative melt 
slope involve either open crystalline structures, such as 
water and graphite, or in the case of closed packed solids, 
a promotion of valence electrons to higher orbitals upon 
compression (6s to 5d in cesiumii, for example). In these 
cases, the liquid is denser than the solid when they coex- 



ist, possibly because of structural or electronic transitions 
taking place continuously in the liquid, as a function of 
pressure, but only at discrete pressure intervals in the 
solid. 

In contrast, recent experiments- have shown that hy- 
drogen phase I - a solid structure with rotationally-free 
molecules associated with the sites of a hexagonal closed 
packed (hep) lattice - persists below the liquid transition 
up to at least 150 GPa, i.e. well beyond the melt curve 
maximum predicted here. The promotion of electrons to 
higher orbitals in hydrogen can also be ruled out because 
of the prohibitive high energy involved. Alternatively, 
it has been suggested that dissociation in the fluid, ei- 
ther gradual or following a first-order liquid-liquid (LL) 
phase transitioni2*i£, may be the origin of a maximum in 
the melt line. This idea is not supported by our results. 
Instead, we explain the physical origin of the maximum 
in terms of changes in the intermolecular interactions - 
a mechanism significantly different from that expected 
from familiar phcnomenological models. 

The approach taken here to compute the melt line is 
one of direct simulation of the melting process. Hystere- 
sis effects of super-heating or super-cooling during the 
phase transition are avoided by simulating solid and liq- 
uid phases in coexistence. The validity of this method 
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FIG. 1: Snapshots from two- phase MD simulations at P = 
130 GPa and temperatures below and above the melting tem- 
perature. A quantitative assessment of the instantaneous lo- 
cal order environment around each molecule has been per- 
formed as described in Ref. ITtI and compared with single- 
phase solid and liquid simulations. Molecules are coloured 
according to the arrangement of their nearest neighbours, red 
and blue representing configurations uniquely characteristic 
of the hep solid and liquid at the given P and T, respec- 
tively (the Q4 order paramete»— has been used for the colour 
map). The systems reach equilibrium at the target average 
temperatures (800 and 900 K) over a time interval of ~ 0.5 
ps, which is followed by periods of coexistence lasting from 1 
to several picoseconds. The phase transitions are observed by 
monitoring changes in the diffusion constants, pair correlation 
functions, specific volumes and local order parameters. In ad- 
dition, some noticeable correlations between temperature and 
diffusion variations indicate the exchange of latent heat. 
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is then assessed by reducing size effects in such a way 
that realistic thermodynamic processes can be mimicked. 
This technique, known as two-phase simulation, has been 
mostly applied with model potentials^, thus allowing 
the simulation of large systems, and only recently it 
was demonstratedi^i^ that implementations within the 
framework of first principles molecular dynamics (MD) 
are feasible. 

Two-phase simulations (Fig. 1) are performed for var- 
ious pressure and temperature conditions to predict the 
melt curve of hydrogen up to 200 GPa (Fig. 2). As dis- 
cussed in Refs. and0, the experimental data available 
up to 44 GPa can be equally well fitted with several em- 
pirical melt equations, leading to qualitatively different 
conclusions regarding the further rise or fall of the phase 
boundary. Our results at 50 GPa compare favourably 
with recent measurement, and the calculations at 130 
and 200 GPa predict a negative slope in this pressure 
region. 

Previous first-principles calculations 12 ' 19 indicated the 
existence of a sharp transition from molecular to non- 
molecular fluid, thus raising the possibility that the melt 
curve maximum is a liquid-liquid-solid triple point. We 
therefore studied the properties of the liquid at temper- 
atures above melting. In our simulations, we do find 



t ■ 1 ■ 1 "t — r 




P (GPa) 

FIG. 2: Melt curve of hydrogen predicted from first princi- 
ples MD. The filled circles are experimental data from Refs. (J 
and and references therein, and the open squares are mea- 
surements from Ref. 0. Triangles indicate two-phase simula- 
tions where solidification (up) or melting (down) have been 
observed, and bracketed melting temperatures (T m ) are repre- 
sented by open circles. As the phase boundary is approached, 
the period of coexistence increases and eventually the out- 
come becomes dependent on the choice of simulation param- 
eters. This degree of arbitrariness is reflected in the error 
bars of T m , which also include the standard deviation of the 
temperatures collected during the MD simulations. All ex- 
perimental and theoretical points are given equal weight and 
fitted with a Kechin melt equation— (solid line in the figure) : 
T m = 14.025 (1 + P/a) b exp(-cP) K, where P is in units of 
GPa, a = 0.030355, b = 0.59991, and c = 0.0072997. The 
open diamond marks the liquid-liquid transition from molec- 
ular to non-molecular fluid at 200 GPa, and the estimated 
slope of this phase boundary is given by the dashed line. The 
error bar on the diamond symbol indicates the hysteresis ef- 
fects during the simulation of the liquid-liquid transition. 



a first-order liquid-liquid phase transition from molec- 
ular to dissociating fluid, but this occurs at tempera- 
tures higher than those of the melting line; at 200 GPa, 
the liquid-liquid transition takes place between 900 and 
1,000 K (Fig. 2). Thus, in the entire pressure range con- 
sidered here, the liquid remains molecular below 900 K, 
with molecules stable over the simulation times of sev- 
eral picoseconds. We emphasize that on heating above 
a critical temperature, the transition to a non-molecular 
fluid takes place very rapidly, over a couple of hundred 
femtoseconds of simulation time. When the fluid is subse- 
quently quenched, the molecules recombine and the sys- 
tem reverts to the molecular state. Therefore, we con- 
clude that the change from a positive to a negative melt 
line slope happens gradually and we emphasize that it 
is not directly related to molecular dissociation. How- 
ever, extrapolations of the melt line and the liquid-liquid 
phase transition indicate a triple point at ~ 300 GPa and 
400 K. Above this pressure, the solid is expected to melt 
into a metallic liquid. 

To estimate the errors originating from finite system 
size effects, we have computed the ground state pressure 
and energy for different size simulation cells sampled at 
the the T-point and with a 2x2x2 k-point mesh. The 
pressure and energy differences between the solid and 
the liquid are converged to within 0.3 GPa and 1 meV 
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FIG. 3: Difference of the specific volumes (AV) and en- 
thalpies {AH) between the liquid and solid phases at the 
melting temperatures determined from the two-phase simu- 
lations. The reported uncertainties include standard devia- 
tions collected during the MD simulations and observed size 
effects from singe-phase simulations with 360 and 768 atom 
supercells. The data are used to compute melt slopes from 
the Clausius-Clapeyron equation: dTm/dP — T m AV/ AH; 
the values obtained for 50, 130 and 200 GPa are (6.5 ± 1.2), 
(-1.4 ±0.6) and (-2.3 ±0.6) K/GPa, respectively. For com- 
parison, the slopes from the melt line fit in Figure 2 are 3.9, 
-2.2 and -2.7 K/GPa, but we note that these are weighted 
heavily by the experimental data, especially at lower pres- 
sure. 
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per ion, and the estimated errors are included in Fig. 3. 
A question that remains to be answered is whether the 
size of the system considered here is sufficient to simulate 
liquid-solid coexistence. To address this, we carried out 
MD simulations of single-phase solid and liquid hydrogen 
at the melting temperatures determined in the two-phase 
simulations. The differences in the specific volumes and 
enthalpies obtained in this way (Fig. 3) are then used to 
compute the melt slopes using the Clausius-Clapeyron 
equation. The agreement with the results from two-phase 
simulations gives us confidence that the computed melt 
curve is thermodynamically consistent. 

We also discuss the validity of the principal approx- 
imations in our theoretical method: (1) adiabatic in- 
teractions, (2) the generalized gradient approximation 
(GGA) for the exchange-correlation energy, and (3) clas- 
sical treatment of ion motion. Electron-phonon inter- 
actions are significant in the solid phase when the di- 
rect band gap is comparable with the phonon frequen- 
cies; an event that eventually occurs at high densities. 
But at 200 GPa and 600 K - the conditions examined 
here where the electron-phonon coupling is expected to 
be most pronounced - the GGA band gap is still ~ 2 eV. 
We note that the GGA tends to systematically underesti- 
mate band gaps; this has been shown for various materi- 
als, including hydrogen 20 . However, because the highest 
vibrational frequencies in hydrogen are ~ 0.5 eV, the 
computed GGA band gap - a lower bound to the ac- 
tual value - is already sufficient to establish that non- 
adiabatic effects are indeed negligible. 

Local density approximations are also expected to 
favour the phase with more delocalised electrons; previ- 
ous studies have shown a consistent tendency of the GGA 
to underestimate melting temperatures, such as those of 
lithium hydride^ and aluminium 16 . Here, this effect is 
expected to be very small because the molecular bond- 
ing properties and the atomic coordination differ little in 
the solid and liquid phases. This also implies that the 
zero-point energy bound in the vibron motion is similar 
in the two phases. Indeed, velocity- velocity autocorrela- 
tion analysis carried out on our MD trajectories shows 
that vibron frequency distributions (containing all vibra- 
tional modes) differ by less than 10 cm -1 in the solid 
and the liquid between 50 and 200 GPa. In addition, wc 
have computed the first-order quantum correction to the 
ionic free energies. It is given by the Wigner-Kirkwood 
formula, AF = h 2 /(24k%T 2 ) J2i( F i)/ m ii where the av- 
erage is over the classical ensemble, and Fj and m, are 
the ionic forces and masses. The difference in the quan- 
tum corrections for the two phases obtained in this way 
at pressures of 50, 130 and 200 GPa remains less than 2 
mcV for temperatures near the computed melt curve. 

We now turn our attention to the physical origin of 
the maximum in the melt curve of hydrogen. The tran- 
sition from an ordered to a disordered state is governed 
by the relative importance of entropy and energy. At 
high pressure, the interactions become mostly repulsive, 
and it is reasonable to expect that a deviation from a 
crystalline arrangement becomes progressively more un- 
favourable in terms of enthalpy. This is indeed the case 
for most materials, as is reflected in the positive slopes of 
their melt curves. On the other hand, attractive many- 
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FIG. 4: The probability distribution of MLWF spreads^ at 
T=700 K and P=50 (black curves), 130 (red curves) and 200 
GPa (green curves). The solid and dashed curves correspond 
to distributions collected from the solid and liquid phases, 
respectively. The inset shows the effective hydrogen-hydrogen 
potentials obtained directly from the simulation trajectories 
with the force-matching procedure described in Ref. (the 
same coloring has been used). As the pressure is increased, 
the spreads of the MLWF are shown to increase faster in the 
liquid than in the solid phase. In addition, effective potentials 
soften more quickly in the liquid than in the solid phase at 
the nearest neighbor molecular separation distance. Note that 
the potentials at 50 GPa in the liquid and solid phases, and 
at 130 GPa in the solid phase, are indistinguishable on this 
scale. 



body interactions can soften the steep increase of the 
repulsive forces at high densitySL 2 ^. Although, in princi- 
ple, such a softening opens up the possibility for a melt 
curve maximum^, it does not appear to be a sufficient 
condition for its existence. Indeed, softening of the re- 
pulsive part of the pair potentials has been observed for 
a number of materials for which there is no experimen- 
tal indication of a melt line turnover. As will be shown 
below, the physics behind the unique melt curve of hy- 
drogen is related to the different rate at which the repul- 
sive interactions (at intermolecular range) are weakened 
in the solid and liquid phases as a function of pressure. 

To investigate how many-body interactions give rise to 
repulsive force softening in the liquid and the solid, we 
have performed an analysis based on maximally-localized 
Wannier functions (MLWF). These functions have been 



shown4» to provide an intuitive description of chemi- 
cal bonds in condensed phases and they have recently 
been used 2 ^ to investigate the infrared activity of candi- 
date structures of hydrogen phase III. We have computed 
MLWF along the 700 K isotherm, at various pressures, 
and find that their spreads (Fig. 4) increase faster in the 
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fluid than in the solid, as the system is compressed. Im- 
portantly, the differences in the two phases come from 
the tails of the MLWF, while their profiles around the 
centres of the molecular bonds remain nearly the same. 
Therefore, the large MLWF spread observed in the high- 
pressure fluid comes from the increased overlap of the 
molecular orbitals, here represented by the MLWF. Be- 
cause the specific volumes of the solid and the liquid 
are very similar under the conditions examined here, the 
physical origin of the larger spread in the fluid can be at- 
tributed to the disorder in the liquid phase. Furthermore, 
the tails of the MLWF in the fluid are characterized by 
an asymmetry, which can be viewed as an effective in- 
termolecular charge transfer, giving rise to an enhanced 
infrared activity in the liquid phase. 

The properties of the MLWF found in our calculations 
are indicative of the behaviour of repulsive interactions in 
the solid and liquid as a function of pressure. Effective 
hydrogen-hydrogen pair potentials derived from our ab 
initio MD trajectories using a force matching technique 26 
are displayed in the inset of Fig. 4. As in the MLWF be- 
haviour, we observe a stronger softening of the potential 
in the liquid than in the solid. Our Wannier function 
analysis indicates that this increased softening is the re- 
sult of charge transfer processes, enhanced by disorder, 
which appear to have a qualitative character similar to 
those proposed in Ref. |2?] for phase III of solid hydrogen. 
We also note that the attenuation of repulsive interac- 
tions observed here in the solid phase correlates well with 
the analysis of Raman spectroscopic data2&, which has 
indicated softening of the effective intermolecular force 
constants above ^100 GPa. 

To summarize, we have predicted the melt curve of 
hydrogen between 50 and 200 GPa from first principles. 
Above approximately 82 GPa, the melt line has a nega- 
tive slope - a phenomenon that we have related to soft- 
ening of the intermolecular interactions, which occur at 
a faster rate in the liquid than in the solid, as a func- 
tion of pressure. Our results provide strong evidence 
that above ^300 GPa the solid melts into a metallic liq- 
uid and that a low-temperature metallic quantum fluid 
will exist at pressures near 400 GPa. Although we have 
only directly probed temperatures as low as 600 K, our 
conclusions are likely to hold at lower temperatures as 
well. Indeed, quantum ion motion usually favours phases 
with lower ion-ion correlations, which would further sta- 
bilize the low-temperature liquid phase found here. In 
addition, our findings are consistent with pressure esti- 
mates from experiments where hydrogen is expected to 
metalliz e - 3 ! 29 . Finally, we note that the observed increase 
in MLWF spreads when going from the high-pressure 
solid to the molecular liquid phase points at the pres- 
ence of dynamically-induced intermolecular charge trans- 
fer (dipole moments); this indicates that the transition 



from solid to liquid at high pressure is accompanied by an 
increase of the hydrogen infrared absorption. Therefore, 
we propose that measurements of infrared activity could 
be used to verify experimentally the predicted turnover 
of the melt curve. 

Methods 
Molecular dynamics 

Our ab initio MD simulations are carried out with the GP 
code, written by F. Gygi, and using the Car-Parrincllo 
(CP) method^ .; the many-body electron problem was 
solved quantum-mechanically within the generalized gra- 
dient approximation 31 of density functional theory, while 
the ions are propagated classically (the validity of this 
approximation for the problem at hand is discussed in 
the text). The actual calculations are performed for deu- 
terium instead of hydrogen. The reason for this is that 
in the classical-ion picture, for the quantities of relevance 
here there is no distinction between hydrogen and its iso- 
topes, but the substitution with heavier ions allows us to 
integrate the equations of motion with reasonably large 
time steps: 2 and 3 a.u. (1 a.u. = 0.0242 fs) depending 
on the density. Furthermore, we use a T-point sampling 
of the Brillouin zone, Troullier-Martins pseudopotentials, 
and 60 Ry energy cut-off for the plane wave expansion of 
the Kohn-Sham orbitals. The accuracy and transferabil- 
ity of our pseudopotential has previously been extensively 
testedi 9 -. We have explicitly verified that the electrons 
remain near the Born-Oppenheimer surface throughout 
the simulations by comparing pressures determined in 
the CP runs with those determined by fully optimizing 
wavefunctions for the same ionic trajectories. The use 
of CP-MD results in a systematic shift of about 0.5 GPa 
both in the solid and liquid phases. 

Two-phase simulations 

We perform the two-phase simulations by equilibrating 
solid (in the hep structure) and liquid phases separately 
in 360-atom supercells with fixed geometry and applied 
periodic boundary conditions. The atomic arrangements 
obtained in this way are then merged, and the liquid 
and solid parts of the resulting 720-atom supercells are 
also heated independently to reduce the strain near the 
newly-formed interface. The whole equilibration process 
typically lasts for several picoseconds of simulation time. 
The MD runs of coexisting phases are then performed 
in the N — P — T ensemble (constant number of parti- 
cles, pressure and temperature respectively). The actual 
system being simulated initially represents semi-infinite 
alternating slabs (by applying periodic boundary condi- 
tions) of solid and liquid, but by the end of the runs only 
one of the two phases, the one that is stable at the chosen 
P and T, fills the entire volume. 
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